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An implementation of both finite rate and equilibrium chemistry have been completed 
for the OVERFLOW code, a chimera capable, complex geometry flow code widely used to 
predict transonic fiowflelds. The implementation builds on the computational efficiency 
and geometric generality of the solver. 


Nomenclature 

y pc 2 /p 

e internal energy (J Q c v dT) 

T| Self Similar Wall Distance for Wedge ((2 — 

z^v)/x) 

Backward Stoichiometric Coefficient of species 
s in reaction r 

Forward Stoichiometric Coefficient of species s 
in reaction r 

0 n Spatial derivative in wall normal direction 

p total mass per unit volume(£™* p s ) 

p s mass of species s per unit volume 

£, Dimensionless Axial Distance for Wedge(x/L) 

c Frozen Sound Speed, dpp^ y 

c s Mass fraction of species s, p s /p 

hr stagnation enthalpy (ei + p/riio + (u 2 + v 2 + 

L Wedge Length(x/L) 


Introduction 

Chemically reacting flows are of increasing impor- 
tance to NASA’s terrestrial and space missions. In 
cases where the flow geometry is complex, which in- 
cludes most practical devices, ability to utilize overset 
grids greatly enhances NASA’s ability to accurately 
model the flow physics. The OVERFLOW code is 
therefore an ideal platform for implementing chemi- 
cal modeling capability. This provides the ability to 
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predict performance and improve design of both reen- 
try vehicles and propulsion systems. In this paper 
we discuss the approach for incorporating finite rate 
chemistry capability into OVERFLOW and present 
representative computations along with comparisons 
with experiment. 

Method 

The base code used is OVERFLOW 1 which is a 
finite difference, chimera(overset) grid capable, com- 
plex geometry flow code widely used for perfect gas 
prediction. Various ADI implicit treatments ( Beam- 
Warming, Pulliam-Chausee diagonal, Symmetric LU 
Gauss Seidell) are available, along with a selection of 
spacial discretization options ( central differnce with 
scalar or matrix dissipation, Roe MUSCL upwinding, 
Yee’s Symmetric TVD, etc.). Multigrid convergence 
acceleration is available, and improves the convergence 
of high Reynolds number viscous flows substantially. 
The dual time scheme implemented in the perfect gas 
version of OVERFLOW allows time accurate calcula- 
tions to be free of factorization and chimera interface 

Tag’. 

Matrix dissipation combined with central space dif- 
ferencing and multigrid was shown to provide a rapid, 
robust solution method suitable for high speed, perfect 
gas flow with strong shock waves. 2 The first step in the 
generalization of the gas model of overflow was to im- 
plement the premixed chemical equilibrium gas model 
of Liu and Vinokur. 3 Preliminary experiments with 
this implementation verified that the central/matrix 
dissipation methodology continued to work with more 
general thermodynamics. 4 

OVERFLOW was modified to replace the perfect 
gas model built into the code with a model assum- 
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ing a gas made up of a mixture of thermally perfect 
(pi = pTIT /M)) gases. This modification affects a 
large fraction of the code, which used perfect gas as- 
sumptions to obtain pressure and temperature from 
the conservative field variables. Four extra field vari- 
ables are added, corresponding to pressure, tempera- 
ture, Prandtl number, and gas constant. 

Chemically Reacting Reynolds Averaged Navier 
Stokes Equations 

The Reynolds averaged Navier-Stokes equations, 
written in conservation law form are 


0Q 9F 9G 9H 

— 4 1 1 

St 3x dy 3z 


- 0 


Where 

Q=[p» pu, pv, pw, pe,] 


F= 


pu, pu z + p + T xx , puv + T X y , puw + T xz 


, , , 0T 

pu(e T + p + T xx ) + VT X y + WT XZ + k— 


dx 
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PV, pVU + Ty X , PV Z + P + Tyy , PVW + Ty 2 , 

, . , 0T1 

pv(ej + P + Tyy ) + WTy z + UTy X + k~ 
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pu, pwu + T zx , pwv + T Z y , pw z + p + x 2z , 


, . , 0T 

pw(e T + p + T zz ) + ut 2X + VT Z y + k— 
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(2) 

(3) 

(4) 

(5) 


The Navier-Stokes equations for a chemically reacting 
gas are similar to this set, except that the 1st equation 
(mass conservation) is replaced by n s equations, one 
for each chemical species. 

0tPi + 0x (upi + (D0 x pi)) + 0y (vpi + (D0yPi)) 

+0 z (wp i + (D0 z p i )) =w i (6) 

These equations add a mass diffusion term 
(0£,D (0£,Pi)) and a source term vOi. The stan- 
dard form for the source terms is 

- v' rt ) tjA. (7) 
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where the r sums the contributions of the individ- 
ual reactions, k f is the forward reaction rate and K* 
is the equilibrium constant(a pure function of tem- 
perature), p r ef is the pressure reference for K e , and 
Av r = T* s w f - Y^ s v b 

r Z_s v rs Z_S rs* 

In this form, the significance of K e is obvious. In 
the case where the mass fractions satisfy the relation 


P 

Pref 




*n(£) =n 


Ps 

Ms 


(9) 


for reaction r, the source term for this particular 
reaction will be zero. K e can be calculated from 
equilibrium thermodynamics, or can be specified in- 
dependently. Both options are available in the current 
implementation . 

In the current implementation, the original Navier- 
Stokes equations(l-S) are solved, with the species frac- 
tions held constant, The species equations (6) are 
then solved holding the velocities and temperatures 
from the mean flow equations just calculated. This 
over specifies the problem, as the mass conservation 
equation is contained in the species mass conservation 
equations. Although various methods of removing this 
over-specification were attempted, the current method 
simply ignores the over-specification, and utilizes the 
Pi to get the partial pressures, and p as the total mass 
density. As the solution converges, the ratio Pi/p 
is used to verify convergence along with the equation 
residual levels. 

Generalization of the Thermodynamic Model 

At the start of each step, just before the viscos- 
ity is calculated, the temperature is found for each 
point from the current internal energy and species 
concentrations, via a Newton-Rap hson iteration. This 
temperature can then be used to find the pressure at 
each point, by summing the partial pressures of each 
gas (p = Lspecies Pspecies). The Prandtl number is 
calculated from Euken’s relation. 5 The gas constant is 
calculated from the species fractions and their molec- 
ular weights. 

The internal energy and absolute entropy of each 
species as a function of temperature are input as a 
table to generalize the thermodynamic model. The 
table has the following format: 

N t To AT 
( M AM 8 1 
M! AH f1 

T 0 t 0 /lU %/K 

Ti fei/ftT h/n 


species 1 


species 


species N t 


y N t 

(mams 2 


Mi 

AHf 2 

To 

t 0 /m %/n 

Ti 

€^/lU f, m 

\Tn t 

s Ut /k 

(MAMS Ns 

•Mn s 

AHf Ns 

To 

e 0 /KT § 0 /n 

Ti 

ti/m $i /n 

VTnt 

tu T /TZJ ?n t /^ 


The first line gives the number of (equis- 
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paced) temperature values used(Nr), the starting 
temperature(T 0 (K°)), and the temperature spacing 
(AT(K°)). The second line is chemical formula for 
the species, such as H20 for water, 02 for diatomic 
oxygen, etc. The next line is the molecular weight 
in kg/kmole and the heat of formation in K Joules /kg. 
The next Nr lines are the values of T(K°), internal en- 
ergy (non dimensionalized by 7£(T), and entropy (non- 
dimensionalized by 7£). This is then followed by a 
similar set of lines for the each species in turn. The 
order of the species in this file is the order of the species 
stored in the solution file, and the order of species for 
specifying reactions rates, etc. 

This choice of functional form for the table, and 
hence the cubic spline fits also results in a smaller 
variation over in the functional values (at least for 
e) over large temperature ranges. A species with 
constant specific heat would have the e/lZl column 
a constant. Tables have been created for 5 species 
air ([N 2> 0 2 ,N0,N,0], from 3 * ), and hydrogen-air 
([H 2 , N 2j 0 2 , HO, H, O, H 2 0, H0 2) H 2 0 2 ], from NIST 
polynomial fits). 

Polynomial fits, such as the Lewis, or NIST curve 
fits may be input by creating a table by evaluating the 
polynomials and writing out the values of T, t/{1t T 
and $/H. Energy and absolute entropy are both input, 
so that chemical equilibrium or finite rate chemistry 
may be calculated. This also allows the backward rate 
for the chemistry to be calculated from the forward 
rate and the 2 nd law of thermodynamics. 

The field equations for convection of species are also 
available in the perfect gas version of OVERFLOW. 
These were replaced with a generalized upwind or 
central convection, 2nd order mass diffusion, and gen- 
eralized chemical source module. These are currently 
solved uncoupled from the fluid and turbulence field 
equations. A fast source calculation procedure was 
implemented, so that the source flux and Jacobian re- 
quires about the same amount of computational effort 
as the efficient P ulliam- Chausee diagonal fluids step. 

Chemical Reaction Specification 

The chemical reactions are input in the OVER- 
FLOW input file, in a new namelist section in the 
input file $CHEM. The number of reactions NRE- 
ACT , the number of elements NELEM, the dis- 
tribution of elements amount the various species 
scoef(ispecies,jelem), the stoichiometric coeficients 
vnuf(ireact,ispecies), xmub( ireact, ispecies, and the re- 
action forward find backward rates fcoef(:, ireact), 
ecoeff:, ireact) are all specified in this section. William 
Chan has written a version of the overgrid 6 which fa- 
cilitates the input of these variables. This software, 
reads the thermodynamics input files gas.tbl described 
above, and parses out the scoe/ matrix from the species 
names, and then allows reactions to be specified in 
an intutive manner. The system also checks for ele- 


ment balance among the chemical reactions input, and 
requires the reactions to balance elements between for- 
ward and backward sides. 


A decision was made to consistently use SI units 
when dimensional values were necessary, such as for 
Too, Poo, reaction rates, etc. The units used then 
are kmoles, m, kg, K° and sec. A great deal of the 
chemical reaction data is in cgs (moles, cm, g, and 
sec) units and factor of 10 -3 will transform a re- 
action rate from cm 3 /(mole sec) to m 3 /(kmole sec). 
The forward rate is input in the Arrenhius form kf = 
Ci (T) C2 cxp(C 3 /T), and the backward rate is deter- 
mined either from an equilibrium coefficient ke = 
kb /kf = k e (ecoefi,T) in the modified Arenhius form 


..nr* J nn 7 r. C . ■ rs 

UOCU 1AX iCiCi — - - 


- +l.n On<4 Inn. 4* V» ni. 
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modynamics. 


The reactions are specified separately, so each that 
third body reactions is input separately. This greatly 
expands the size of this input section, but allows for 
complete generality in fixing the forward and backward 
rate of each reaction for the various third bodies. The 
chemical source calculation procedure goes through 
the complete reaction set, and groups reactions from 
the same third body reaction together, taking account 
of the variations in the forward and backward rates 
among the group of various third body species, for 
computational efficiency. 


For the species equations, both central and 2nd or- 
der upwind are available, and solutions using both 
approaches are given in the results section. There are 
cases for which the upwind method is more robust, but 
once initial transients are eliminated, the central space 
version is more robust ( can run at higher CFL) and 
accurate ( as shown by the element conservation). 


Transport Properties 

The diffusion operator for the species equations is 
2nd order accurate, and the binary diffusion coeffi- 
cients and laminar viscosities are given by Leonard- 
Jones theory, with a Wilke mixing rule. Euken’s rela- 
tion is used to give the molecular thermal conductivity, 
or Prandtl number, of the gas mixture. Turbulence is 
modelled by any of the available turbulence models 
in OVERFLOW. The effect of turbulence is only to 
increase the effective viscosity, binary diffusivity, and 
heat transfer. No enhancement of the reaction rates is 
currently modelled. 


Computational Efficiency 

A finite rate simulation is a little more than 3 times 
as expensive per iteration as a corresponding simu- 
lation using the perfect gas version of OVERFLOW. 
The additional work ( with a perfect gas simulation 
as a unit of computational work), is as follows: 
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Arbitrary Thermodynamics (single species) 0.25 

Species Convection, etc. (5 species) 1.5 

Source Flux/ Jacobian (5 species) 0.5 

Note that the species convection/diffusion includes the 
pre-multiplication by a 5 x 5 matrix, required for the 
ADI solution of the species equations with implicit 
source terms. This is substantial work compared to 
the ADI inversion itself, and makes the species field 
operator an expensive part of this good breakfast, er. 
flow solver. 


Non-Premixed Chemical Equilibrium Capability 

Along with the new finite rate capability, an ability 
to accurately simulate flows in chemical equilibrium 
has been implemented. This capability was envisioned 
as a means of providing a robust starting solution ca- 
pability for difficult flowfields. This capability does not 
require that the flow be premixed, so that equilibrium 
diffusion flame calculations are possible. For flowfields 
where chemical equilibrium is a valid assumption, this 
provides a robust and more simply specified simulation 
capability. 

This capability is provided by a similarity transfor- 
mation of the species variables. The species can be 
thought of as composed of linear combinations of the 
elements. For instance, water (H 2 0) can be thought 
of as 2 parts hydrogen to one part oxygen (on a molar 
basis). A linear combination of the species then can be 
used to compute the molar concentraion of each ele- 
ment, given the species concentrations. For a chemical 
system composed of rt e elements, and n s species, there 
are n e linear equations which will give the amount 
of each element given the species concentraions. Fur- 
thermore, on physical grounds, these rt e equations are 
independent, and thus define a n e dimensional sub- 
space within the n s dimensional species space. A 
method equivalent to Gram-Schmidt orthogonaliza- 
tion can then generate the remaining n s — n e linear 
relations which will transform the original species vari- 
ables into a new combinations where rt e of the new 
variables correspond with the elements, and the re- 
maining n s — rt e variables correspond to the degrees 
of reaction freedom of the system. 

This transformation can (and has been) automated, 
and has some desirable properties which merit explo- 
ration. One of these properties is that by transforming 
the species equations to this new independent variable 
set separates out the rt e zero eigenvalues of the source 
jacobian matrix. More than this, these variables have 
zero source, so they could be separated from the other 
rt s — Tie variables, and solved as simple convection 
diffusion equations. The other ’reaction freedom’ vari- 
ables can then be solved as the original species set, 
with non-zero sources, and a matrix block size reduced 
by the n e elements sepatated out at considerable com- 
puational savings. 

Currently, this ’’element splitting” is utilized to im- 


plement the non-premixed chemical equilibrium ca- 
pability of the code. The species are split, and the 
equations for the element conservation are solved. The 
’reaction freedom’ variables are solved for by enforcing 
chemical equilibrium at each point. 


Multigrid, Relaxation, and Chemistry-Fluid 
Coupling 


The multigrid method 8 has been modified to en- 
hance the robustness of the method. The original 
multigrid method used a linear interpolation of the 
conservative variables to map the coarse grid solu- 
tion onto the next finer grid. This interpolation was 
replaced with a linear interpolation of the variables 
[p,u,v, w, e], which ensures a positive value of e on 
the finer grid if e is positive on the coarser grid. With 
this modification, the multigrid method continues to 
be a net win, enhancing robustness and improving con- 
vergence rates. 


The Pulliam-Chaussee diagonal algorithm is used 
for the mass, momentum and energy equations, and 
enhances the computational efficiency. This implicit 
operator requires that the eigensystem (here essen- 
tially just the sound speed) be modelled correctly. The 


frozen sound speed 
with 7 == 1 + <*7 • 


dp 


S,Ci 


is calculated by (yp/p), 


The chemistry implementation is currently a loosely 
coupled approach, with the changes in mass, momen- 
tum and energy variables {p, pu, pv, pw, pe} solved for 
holding species concentration fixed, and then changes 
in species concentration {Ci} are determined by solving 
the conservative form of the species transport equa- 
tions holding the conservative variables constant. 


Boundary Conditions 

A new boundary condition was implemented for 
chemically reacting flowfields, the catalytic fixed tem- 
perature wall. The non-catalytic boundary condition 
for the species is zero gradient in species fractions. A 
catalytic wall boundary condition is given as 


WMpn,)=PnYc 

WWpnI—PnT. 

ra „( P 0 ! )= P 0 T .(^) 1 

TO.(po)=-pot c (^) T 

£>9n(p No)=0 


( 10 ) 


where the value of y c is 0 for a non-catalytic wall, or 
1 for a fully catalytic wall. 

This is implemented in a general way by defining 
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the diffusion potential of the species j as 





1 

2 


then we can write a general version of this boundary 
condition as 


2?0n( P j) = CjkWlc (11) 

k 

where Cjk is a n s x ns matrix. 

Results 



Fig. 1 10°M=25 Inviscid Wedge Grid, Colored by 

Axial Velocity 


Inviscid Wedge Flow 

This flowfield is a self similar flowfield useful for 
gauging the ability to capture shocks at oblique an- 
gles to the grid, and to obtain the external inviscid 
flowfield downstream of the shock. The case consid- 
ered here is a 10° wedge, at M=25.083, with freestream 
static conditions of p**, = lOOKPa,^ = 2000° K. The 
species for this case are {N 2 , O 2 , NO, N, O}. The ther- 
modynamic data for these species (cv(T)) are from the 
calculations outlined in Liu. 3 

A converged solution is shown in Figure 1. The cap- 
ture of the oblique shock is clearly demonstrated, with 
uniform conditions obtained on both sides of the shock. 
The shock crosses the grid system obliquely, with the 
shock moving out one grid line for every two gridlines 
moved downstream. 



Fig. 2 Inviscid M=25 Wedge, Profiles of Pressure, 
Total Enthalpy, Temperature and y 



Fig. 3 Inviscid M=25 Wedge, Mass Fraction Pro- 
files 



Fig. 6 Circular Cylinder Flowfield: Moo = 12.7, 
r = 1m, h — 40km, |v| = 4.04km/s 


Profiles of pressure, total enthalpy, temperature, 
and pressure are shown in Figure 2. The pressure 
profile shows minimal overshoot, with good similarity 
between the various axial locations The total enthalpy 
is well conserved (for this flowfield it is constant), with 
errors, even at the shock, less than 0.5%. Total en- 
thalpy errors behind the shock are less than 0.1%. 
The post shock temperature reached in this flowfield is 
8000° K, for which there is appreciable dissociation, as 
shown by the gamma variation from 1.32 (freestream) 
to 1.26(post shock). 

The species mass fractions are shown in Figure 3. 
Almost all the diatomic oxygen (O 2 ) is dissociated 
across the shock, and is transformed into nitrous oxide 
(NO). Elemental oxygen (O) and nitrogen (N) are at 
low levels post shock. 

The forward and backward rate for this calculation 
are given by Parks model, 7 and in the full paper this 
will be compared with solutions using backward rates 
derived from the forward rates and the equilibrium 
constants given by the 2 nd law. 

2D Cylinder Flow 

As a comparison to other Codes, a solution for a 
lm radius circular cylinder at 40km altitude, Uoo = 
6 km/sec (figure ) is compared with solutions from 
GASP. The flow conditions for this case are = 
12.7, Re r = 9690, = 251.05°K, and V co = 277.5Pa. 

The freestream mass fractions are cn 2 = 0.767, Co 2 = 
0.233 with traces of NO, N, and O. Two wall boundary 
conditions are evaluated, non catalytic and catalytic 
fixed temperature, with T w = 2553° K. 

The dimensions of the finest grid were 101 x 101 , and 
the solution was evaluted on coarser grids of 51 x 51 
and 26 x 26 to assess grid convergence. The coarser 
grids were obtained from the finer grid by successively 
deleting every other point. Wall normal spacing for the 
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finest grids was 3x10 5 for both GASP and OVER- 
FLOW. 

One reaction model for this case was created by us- 
ing the forward rates of the Park model, and backward 
rates from equilibrium thermodynamics. A second 
model uses Park’s equilibrium coeffiecient curve fits 
for n = le!8, as K e (T), then setting p re f = lOOPa. 
These solutions are compared with the solution from 
GASP, using its coded Park model. 

Profiles of the thermodynamic variables along the 
stagnation streamline are shown in figure 5. Total 
enthalpy is well conserved across the shock, and the 
predictions of GASP and the two reaction models 
are in good agreement. Pressure density are also in 
very good agreement, and the OVERFLOW solution 
appears to be grid resolved on the 51 x 51 grid. Tem- 
perature however, is not in as good agreement as the 
other variables. 

The reason for this can be seen from figure 4, where 
the stagnation line species concentrations are shown. 
The levels of monatomic nitrogen and oxygen are well 
predicted, but there is a slight disagreement predic- 
tions for NO, N 2 , and 0 2 . These can be traced to 
differences in the K e representations, as the results can 
be brought into much closer agreement by altering the 
reference pressure for K e . As there is some slop in the 
various estimates of K e for these reactions, this agree- 
ment is probably as good as can be expected. 

F4 Nozzle Flow 



Fig. 7 T10-97 Pressure Predictions, Finite Rate 

Another important test case for chemically react- 
ing flows are nozzle flowfields. 'Test 10-97(figure7) of 
the HHSFD database is the flow in the ONERA F4 
wind tunnel nozzle. This flowfield utilizes synthetic 
air, containing only nitrogen and oxygen ( no argon, 
C0 2 , etc.), with mass fractions of 7/9 and 2/9 respec- 
tively. The stagnation chamber of the flowfield has 
extremely high enthalpy (100 times standard sea level 
conditions) and extremely high pressure(37.3 MPa), 
The flow then expands out a axisymmetric nozzle to 
near vacuum conditions (p=52Pa). The flowfield has 
been computed with a variety of codes, 9 providing a 
useful comparison with other codes, as well as experi- 
mental data. 

Inlet conditions are highly dissociated subsonic flow. 
Plug flow in chemical equilibrium was imposed at the 
inlet face of the flowfield, and p t and T t were cho- 
sen to produce the correct enthalpy and pressure in 
the stagnation chamber. The Spalart-Almaras turbu- 
lence model was used, and the initial -vj was set to 
0.1 laminar viscosity. Tunnel walls are modelled as 
constant temperature no-slip boundaries, and the exit 


plane boundary condition is simple extrapolation. 

The flow is close to being in chemical equilibrium, 
and equilibrium codes are as successful as finite rate 
codes in predicting the nozzle wall pressures. Figure 
8 shows a comparison of the wall pressure predictions 
of the code with experimental values. Also shown are 
the wall pressure predictions of the code in equilib- 
rium mode (using the Parks equil. coefficients), and 
the premixed chemical equilibrium version of OVER- 
FLOW. 4 

Exit pressures are reasonably well predicted for this 
case, as shown in figure 9, and are very much in 


Wall preessures, T10— 97 of RTO Test Cases 


\ 

\T\ 

\\ 

\\ 

\ 

> 

Leger 

o 

(X 
\ \ 

\ \ 
\ ’ 

’ 512x129 
512x129, 
EXP 

jL 

Chemical i 

quffibrium 

- 


■ 

T — 7 

/ / 

■sK 


kf 
. . 


0.5 1 1.5 2 2.5 3 3.5 


x(m) 


Fig. 8 Wall Pressure Field Predictions, Compari- 
son with Experiment 



Fig. 9 Exit Plane Pitot Pressure Field Predictions, 
Comparison with Experiment 
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line with predictions of other nonequilibrium codes for 
this case. The exit temperature is also well in line 
with other predictions, and agrees well within the ex- 
perimental uncertainty (figure 10). A final check is 
provided by the NO mass fraction prediction, as seen 
in figure 11. This prediction agrees well with other 
predictions for this quantity, in terms of the the exit 
level, and how quickly it approaches this value. 

Other Results 

Viscous wall solutions for other geometries includ- 
ing 2D and 3D reentry configurations and nozzles, will 
be presented in the full paper, along with comparison 
with experimental results. 



Fig. 10 Centerline Temperature Predictions 
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Fig. 11 Centerline NO Mass Fraction Predictions 
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Conclusions 

Finite rate and equilibrium chemistry capability has 
been added to a version of the OVERFLOW code. The 
implementation is capable of predicting both combus- 
tion and reentry flowfields, and extends the overset 
capability of OVERFLOW into reacting flowfields. 
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